# Replication file for:

# The Effect of Wartime Legacies on Electoral Mobilization after Civil War
# Felix Haass, Martin Ottmann

# Datestamp: 29.09.2021

# Description: This file replicates all tables and figures in the appendix.

# Usage: please follow the instructions in the file README.txt


# Check if packages are installed -----------------------------------------


# Check if packages are installed -----------------------------------------

list.of.packages <- c("here",
                      "lfe",
                      "fastDummies",
                      "tidyverse",
                      "patchwork",
                      "fixest",
                      "lmtest", 
                      "multiwayvcov",
                      "here",
                      "RcppArmadillo",
                      "Rcpp",
                      "geosphere",
                      "ggsn",
                      "stargazer",
                      "kableExtra",
                      "sf",
                      "data.table",
                      "ggpubr",
                      "modelsummary")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)


# Load Libraries ---------------------------------------------------------------

library(here)
library(lfe)
library(fastDummies)
library(tidyverse)
library(patchwork)
library(fixest)
library(here)
library(lmtest)
library(multiwayvcov)
library(stargazer)
library(kableExtra)
library(sf)
library(ggpubr)
library(ggsn)
library(modelsummary)
library(RcppArmadillo)
library(Rcpp)
library(geosphere)
library(data.table)


# Load data ---------------------------------------------------------------

model_data_kec_month <- readRDS(here("./data/model_data_kec_month.rds"))
model_data_kab_month <- readRDS(here("./data/model_data_kab_month.rds"))

vh_clean <- readRDS(here("./data/vh_clean.rds"))
election_2006 <- readRDS(here("./data/election_2006.rds"))


# Summary statistics - RAN Data -------------------------------------------

subdist_sumstats <- model_data_kec_month %>% 
  ungroup() %>% 
  select(
         pre_election_2006, 
         control_dummy_median_nvms, 
         control_nvms, 
         number_of_new_projects, 
         committed_usd, 
         contains("sector"), 
         contains("partner")) %>% 
  select(-sector_multi_count) 

subdist_sumstats <- subdist_sumstats %>% 
  rename(`Aid: Number of New Projects` = number_of_new_projects, 
         `Aid: Committed USD` = committed_usd, 
         `Aid Sector: Administrative` = sector_admin_count, 
         `Aid Sector: Education` =  sector_educ_count, 
         `Aid Sector: Economic Development` = sector_econ_count, 
         `Aid Sector: Health`= sector_health_count, 
         `Aid Sector: Infrastructure` = sector_infra_count, 
         `Aid Sector: Institutional Development` = sector_institution_count,
         `Aid Sector: Religion` = sector_religion_count, 
         `Aid Sector: Social` = sector_social_count, 
         `Aid Sector: Unallocated` = sector_unall_count, 
         `Aid Sector: Spatial Planning and Development` =  sector_spatial_count, 
         `Aid Partner: Private` = partner_private_count,
         `Aid Partner: International private` = partner_international_private_count,
         `Aid Partner: Government` = partner_gov_count,
         `Aid Partner: UN` = partner_io_count,
         `Aid Partner: National NGO` = partner_natl_ngo_count,
         `Aid Partner: International NGO` = partner_intl_ngo_count,
         `Aid Partner: Bilateral` = partner_bilateral_count,
         `GAM support (dummy)` = control_dummy_median_nvms,
         `GAM support (continuous)` = control_nvms, 
         `Pre-Eelection 2006 Dummy` = pre_election_2006)

stargazer(subdist_sumstats %>% as.data.frame(), digits = 2, type = "html", 
          summary = T,
          summary.stat = c("n", "mean", "median", "sd", "min", "max"),
          title = c("Obs", "Mean", "Median", "SD", "Min", "Max"),
          out = "./tables/appendix_table_1_sumstats.html", 
          float = F)


# Summary statistics - Kabupaten data -------------------------------------

# Summary statistics - RAN Data -------------------------------------------

dist_sumstats <- model_data_kab_month %>% 
  ungroup() %>% 
  select(
    pre_election_2006, 
    control_dummy_median_nvms, 
    control_nvms, 
    number_of_new_projects, 
    committed_usd) 

dist_sumstats <- dist_sumstats %>% 
  rename(`Aid: Number of New Projects` = number_of_new_projects, 
         `Aid: Committed USD` = committed_usd,
         `GAM support (dummy)` = control_dummy_median_nvms,
         `GAM support (continuous)` = control_nvms, 
         `Pre-Eelection 2006 Dummy` = pre_election_2006)

stargazer(dist_sumstats %>% as.data.frame(), digits = 2, type = "html", 
          summary = T,
          summary.stat = c("n", "mean", "median", "sd", "min", "max"),
          title = c("Obs", "Mean", "Median", "SD", "Min", "Max"),
          out = "./tables/appendix_table_2_district_sumstats.html", 
          float = F)


# Summary statistics - Survey ---------------------------------------------

survey_sumstats <- vh_clean %>% 
  ungroup() %>% 
  select(election_tsunami_06, 
         election_tsunami_05,
         H136_support0105,
         H060i_grp_kpa ,
         tsunami_affected,
         H027_vill_pop,
         H032a_perc_aceh,
         H036_idp,
         H048_poorHH,
         H053a_primary98,
         number_of_vh,
         H110_conf_victims0105,
         H147_lang,
         Hq152_r_6_KPA) %>%
  
  mutate(H136_support0105 = H136_support0105 == 1, 
         H027_vill_pop = log(H027_vill_pop),
         H032a_perc_aceh = log(H032a_perc_aceh),
         H036_idp = factor(H036_idp), 
         H110_conf_victims0105 = H110_conf_victims0105 > 0, 
         H147_lang_aceh = H147_lang == 1, 
         H147_lang_other = H147_lang == 50, 
         Hq152_r_6_KPA) %>% 
  select(-H147_lang) %>% 
  rename("Majority in village supported GAM" = H136_support0105,
         "Did the village receive tsunami aid in 2006?" = election_tsunami_06,
         "Did the village receive tsunami aid in 2005?" = election_tsunami_05,
          "KPA active in village? (Q60)" = H060i_grp_kpa , 
          "Tsunami affected" = tsunami_affected, 
          "Village Pop. (log) (Q29)" = H027_vill_pop,
          "Percent Ethnic Acehnese (log) (Q32)" = H032a_perc_aceh, 
         "IDPs (Q36)" = H036_idp, 
         "Number of Poor Households (Q48)" =  H048_poorHH, 
         "Primary Schools (Q53)" = H053a_primary98, 
         "Number of Village Heads" = number_of_vh, 
          "Conflict Victims (Q110)" = H110_conf_victims0105, 
         "Interviewer Language - Acehnese" = H147_lang_aceh, 
         "Interviewer Language - Other" = H147_lang_other,
         "Was KPA present during Interview?" = Hq152_r_6_KPA)
         
stargazer(survey_sumstats %>% as.data.frame(), digits = 2, type = "html", 
          summary = T,
          summary.stat = c("n", "mean", "median", "sd", "min", "max"),
          title = c("Obs", "Mean", "Median", "SD", "Min", "Max"),
          out = "./tables/appendix_table_3_survey_sumstats.html", 
          float = F)


# Different NVMS cutoffs --------------------------------------------------

## Figure ------------------------------------------------------------------

level_of_control <- ggplot(model_data_kec_month, aes(x = control_nvms)) + 
 
  annotate("rect", xmin = -0.5, xmax = 0.5, ymin = -Inf, ymax = Inf,
            fill = "#f7f7f7", color = "black") +
  annotate("rect", xmin = 0.5, xmax = 4.5, ymin = -Inf, ymax = Inf,
            fill = "#cccccc", color = "black") +
  annotate("rect", xmin = 4.5, xmax = 9.5, ymin = -Inf, ymax = Inf,
            fill = "#969696", color = "black") +
  annotate("rect", xmin = 9.5, xmax = Inf, ymin = -Inf, ymax = Inf,
            fill = "#525252", color = "black") +
  geom_histogram(fill = "white", color = "black", 
                 binwidth = 1, size = 0.6) +
  geom_vline(aes(xintercept = median(control_nvms)), 
             color = "firebrick", size = 2) + 
  geom_label(data = tibble(x = c(3.5),
                           y = 1450,
                           label = c("Median")),
             aes(x = x, y = y, label = label),
             color = "firebrick",
             inherit.aes = FALSE) +
  scale_y_continuous(limits = c(0, 2000)) +
  geom_label(data = tibble(x = c(0, 2.5, 7.5, 15),
                           y = 1750,
                           label = c("No", "Low", "Medium", "High")),
             aes(x = x, y = y, label = label),
             inherit.aes = FALSE) +
  labs(x = "GAM Support Measure (NVMS)", y = "Count of Observations\nin Kecamatan Data") +
  theme_bw() +
  theme(panel.grid = element_blank())

ggsave(plot = level_of_control, filename = "./figures/appendix_figure_4_control_cutoffs.pdf", 
       width = 10, height = 4)


## Models ------------------------------------------------------------------

# kecamatan
aid_model_kec_categorical <- felm(log(number_of_new_projects + 1) ~ 
                                pre_election_2006*control_categories
                              |  kecid + yearmonth_fct | 0 | kecid , 
                              data = model_data_kec_month   )

# kabubaten
aid_model_kab_categorical <- felm(log(number_of_new_projects + 1) ~ 
                                pre_election_2006*control_categories
                              |  KABNO + yearmonth | 0 | KABNO , 
                              data = model_data_kab_month)

# USD as DV
# kecamatan
aid_model_kec_usd_categories <- felm(log(committed_usd + 1) ~
                                       pre_election_2006*control_categories
                                     
                                     |  kecid +yearmonth_fct | 0 | kecid,
                                     data = model_data_kec_month )

# kabupaten
aid_model_kab_usd_categories <- felm(log(committed_usd + 1) ~ 
                                       pre_election_2006*control_categories
                                     |  KABNO + yearmonth | 0 | KABNO , 
                                     data = model_data_kab_month)

stargazer::stargazer(aid_model_kec_categorical, 
                     aid_model_kab_categorical,
                     aid_model_kec_usd_categories, 
                     aid_model_kab_usd_categories,
                     out = "./tables/appendix_table_5_categories.html", 
                     covariate.labels = c("GAM Area (Low Support) X Pre-Election", 
                                          "GAM Area (Medium Support) X Pre-Election", 
                                          "GAM Area (High Support) X Pre-Election") ,
                     dep.var.caption = "", 
                     dep.var.labels = c("No. of New Projects (log +1)",
                                        "Aid Volume (log + 1)"),
                     omit = c(1:4),
                     add.lines = list( c("Geographical Unit", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{District}", 
                                         "\\multicolumn{1}{c}{Sub-District}"),
                                       c("Unit FE", rep("\\multicolumn{1}{c}{Yes}", 8)), 
                                       c("Month FE", rep("\\multicolumn{1}{c}{Yes}", 8))),
                     float = F, 
                     digits = 2,
                     notes = "",
                     notes.append = F, 
                     notes.label = "", 
                     align = T, 
                     keep.stat = c("n", "adj.rsq"))


# UCDP replication --------------------------------------------------------

kec_model_ucdp <-  felm(log(number_of_new_projects + 1) ~ 
                                  pre_election_2006:control_dummy_median_ucdp 
                                |  kecid + yearmonth_fct | 0 | kecid , 
                                data = model_data_kec_month )

kec_model_usd_ucdp <-  felm(log(committed_usd + 1) ~ 
                          pre_election_2006:control_dummy_median_ucdp 
                        |  kecid + yearmonth_fct | 0 | kecid , 
                        data = model_data_kec_month )

kab_model_ucdp <- felm(log(number_of_new_projects + 1) ~ 
                            pre_election_2006:control_dummy_median_ucdp
                          |  KABNO + yearmonth | 0 | KABNO , 
                          data = model_data_kab_month)

kab_model_usd_ucdp <- felm(log(committed_usd + 1) ~ 
                         pre_election_2006:control_dummy_median_ucdp
                       |  KABNO + yearmonth | 0 | KABNO , 
                       data = model_data_kab_month)

kec_model_ucdp_lpm <-  felm(I(number_of_new_projects > 0) ~ 
                          pre_election_2006:control_dummy_median_ucdp 
                        |  kecid + yearmonth_fct | 0 | kecid , 
                        data = model_data_kec_month )

stargazer::stargazer(kec_model_ucdp, 
                     kab_model_ucdp,
                     kec_model_usd_ucdp, 
                     kab_model_usd_ucdp, 
                     kec_model_ucdp_lpm,
                     out = "./tables/appendix_table_7_ucdp.html", 
                     covariate.labels = c("GAM Area (UCDP) X Pre-Election") ,
                     dep.var.caption = "", 
                     dep.var.labels = c("No. of New Projects (log +1)",
                                        "Aid Volume (log + 1)",
                                        "Pr(New Aid Project > 0)"),
                     omit = "as.numeric|Constant",
                     add.lines = list( c("Geographical Unit", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{District}", 
                                         "\\multicolumn{1}{c}{Sub-District}"),
                                       c("Unit FE", rep("\\multicolumn{1}{c}{Yes}", 8)), 
                                       c("Month FE", rep("\\multicolumn{1}{c}{Yes}", 8))),
                     float = F, 
                     digits = 2,
                     notes = "",
                     notes.append = F, 
                     notes.label = "", 
                     align = T, 
                     keep.stat = c("n", "adj.rsq"))


# Region-specific analysis ------------------------------------------------

west_aceh <- c(1171,  
               1108, 
               1116, 
               1107, 
               1114,
               1112, 
               1103 )

model_west_aceh <- felm(log(number_of_new_projects + 1) ~ 
                          pre_election_2006*aceh_ethnic_settle*control_dummy_median_nvms
                        |  kecid + yearmonth_fct | 0 | kecid , 
                        data = model_data_kec_month %>% 
                          filter(!(KABNO %in% west_aceh)))

model_no_aid_removed <- felm(log(number_of_new_projects + 1) ~ 
                               pre_election_2006:control_dummy_median_nvms 
                             |  kecid + yearmonth_fct | 0 | kecid , 
                             data = model_data_kec_month %>% 
                               group_by(kecid) %>% 
                               filter(sum(number_of_new_projects) > 0) )


# Alternative measures of GAM support in ARLS ------------------------------

mod_support_gov <- felm(election_tsunami_06 ~
                          gam_support_gov_0105 +
                          # structure
                          H060i_grp_kpa  +
                          # demographics
                          I(tsunami_affected > 0) +
                          log(H027_vill_pop) +
                          log(H032a_perc_aceh+1) +
                          factor(H036_idp) +
                          H048_poorHH +
                          H053a_primary98 +
                          number_of_vh +
                          I(H110_conf_victims0105 > 0) +
                          # interviewer FE
                          factor(H147_lang) +
                          factor(Hq152_r_6_KPA)
                        | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

mod_support_forest <- felm(election_tsunami_06 ~
                             gam_support_forest_0105 +
                             # structure
                             H060i_grp_kpa  +
                             # demographics
                             I(tsunami_affected > 0) +
                             log(H027_vill_pop) +
                             log(H032a_perc_aceh+1) +
                             factor(H036_idp) +
                             H048_poorHH +
                             H053a_primary98 +
                             number_of_vh +
                             I(H110_conf_victims0105 > 0) +
                             # interviewer FE
                             factor(H147_lang) +
                             factor(Hq152_r_6_KPA)
                           | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

mod_support_food <- felm(election_tsunami_06 ~
                           gam_food_0105  +
                           # structure
                           H060i_grp_kpa  +
                           # demographics
                           I(tsunami_affected > 0) +
                           log(H027_vill_pop) +
                           log(H032a_perc_aceh+1) +
                           factor(H036_idp) +
                           H048_poorHH +
                           H053a_primary98 +
                           number_of_vh +
                           I(H110_conf_victims0105 > 0) +
                           # interviewer FE
                           factor(H147_lang) +
                           factor(Hq152_r_6_KPA)
                         | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

mod_support_sleep <- felm(election_tsunami_06 ~
                            
                            I(gam_sleep_0105 == "Sometimes" | 
                                gam_sleep_0105 == "Often" | 
                                gam_sleep_0105 == "Always" ) +
                            # structure
                            H060i_grp_kpa  +
                            # demographics
                            I(tsunami_affected > 0) +
                            log(H027_vill_pop) +
                            log(H032a_perc_aceh+1) +
                            factor(H036_idp) +
                            H048_poorHH +
                            H053a_primary98 +
                            number_of_vh +
                            I(H110_conf_victims0105 > 0) +
                            # interviewer FE
                            factor(H147_lang) +
                            factor(Hq152_r_6_KPA)
                          | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

stargazer(mod_support_gov, 
          mod_support_forest,
          mod_support_food, 
          mod_support_sleep,
          omit = "H0|vh|tsunami|lang|Hq|Constant|H110",
          covariate.labels = c("Did Gov't consider this village GAM? (Q 126)",
                               "Was a GAM base nearby this village? (Q 127)", 
                               "Did village voluntarily provide food for GAM? (Q 132)", 
                               "Did GAM soldiers sleep in the village? (Q 128)"),
          out = "./tables/appendix_table_6_alt_survey_measures.html",
          type = "html", 
          dep.var.caption = "",
          dep.var.labels = c("Did village receive tsunami aid in 2006?"),
          add.lines = list(c("Kab. FE", rep("\\multicolumn{1}{c}{Yes}", 4)), 
                           c("Covariates", rep(c("\\multicolumn{1}{c}{Yes}"), 4))), 
          float = F, 
          digits = 2,
          notes = "",
          notes.append = F, 
          notes.label = "", 
          align = T, 
          model.names = F,
          column.separate = c(2, 2, 2, 2), 
          keep.stat = c("n", "adj.rsq", "aic"))


# Unit-specific trends ----------------------------------------------------

# kecamatan
mod_kec_trend <- felm(log(number_of_new_projects + 1) ~ 
                        pre_election_2006:control_dummy_median_nvms 
                      |  kecid_fct + yearmonth_fct + kecid_fct:yearmonth | 0 | kecid , 
                      data = model_data_kec_month %>% 
                        ungroup() %>% 
                        mutate(kecid_fct = as.factor(kecid)))

# kecamatan: committed_usd
mod_kec_trend_usd <- felm(log(committed_usd + 1) ~ 
                        pre_election_2006:control_dummy_median_nvms 
                      |  kecid_fct + yearmonth_fct + kecid_fct:yearmonth | 0 | kecid , 
                      data = model_data_kec_month %>% 
                        ungroup() %>% 
                        mutate(kecid_fct = as.factor(kecid)))


# Weighted by unit size  --------------------------------------------------

# kecamatan
mod_kec_area <- felm(log(number_of_new_projects + 1) ~ 
                             pre_election_2006:control_dummy_median_nvms 
                           |  kecid_fct + yearmonth_fct  | 0 | kecid , 
                           data = model_data_kec_month %>% 
                             ungroup() %>% 
                             mutate(kecid_fct = as.factor(kecid)), 
                           weights = log(model_data_kec_month$area))

mod_kec_area_usd <- felm(log(committed_usd + 1) ~ 
                       pre_election_2006:control_dummy_median_nvms 
                     |  kecid_fct + yearmonth_fct  | 0 | kecid , 
                     data = model_data_kec_month %>% 
                       ungroup() %>% 
                       mutate(kecid_fct = as.factor(kecid)), 
                     weights = log(model_data_kec_month$area))


# Unit-specific trends + Weight by size -----------------------------------

# kecamatan
mod_kec_trend_area <- felm(log(number_of_new_projects + 1) ~ 
                        pre_election_2006:control_dummy_median_nvms 
                      |  kecid_fct + yearmonth_fct + kecid_fct:yearmonth | 0 | kecid , 
                      data = model_data_kec_month %>% 
                        ungroup() %>% 
                        mutate(kecid_fct = as.factor(kecid)), 
                      weights = log(model_data_kec_month$area))

# kecamatan: committed_usd
mod_kec_trend_usd_area <- felm(log(committed_usd + 1) ~ 
                            pre_election_2006:control_dummy_median_nvms 
                          |  kecid_fct + yearmonth_fct + kecid_fct:yearmonth | 0 | kecid , 
                          data = model_data_kec_month %>% 
                            ungroup() %>% 
                            mutate(kecid_fct = as.factor(kecid)), 
                          weights = log(model_data_kec_month$area))

stargazer::stargazer(mod_kec_trend, 
                     mod_kec_area, 
                     mod_kec_trend_area, 
                     mod_kec_trend_usd, 
                     mod_kec_area_usd, 
                     mod_kec_trend_usd_area,
                     type = "latex",
                     out = "./tables/appendix_table_8_trends.html", 
                     covariate.labels = c("GAM Area  X Pre-Election") ,
                     dep.var.caption = "", 
                     dep.var.labels = c("No. of New Projects (log +1)",
                                        "Aid Volume (log + 1)"),
                     add.lines = list( c("Geographical Unit", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{Sub-District}", 
                                         "\\multicolumn{1}{c}{Sub-District}"),
                                       c("Unit FE", rep("\\multicolumn{1}{c}{Yes}", 8)), 
                                       c("Month FE", rep("\\multicolumn{1}{c}{Yes}", 8)), 
                                       c("Unit-specific Trends", rep(c("\\multicolumn{1}{c}{Yes}", "\\multicolumn{1}{c}{No}", "\\multicolumn{1}{c}{Yes}"), 2)),
                                       c("Weights", rep(c("\\multicolumn{1}{c}{No}", "\\multicolumn{1}{c}{Yes}", "\\multicolumn{1}{c}{Yes}"), 2))),
                     float = F, 
                     digits = 2,
                     notes = "",
                     notes.append = F, 
                     notes.label = "", 
                     align = T, 
                     keep.stat = c("n", "adj.rsq"))


# Plot Trend lines -------------------------------------------------------------

did_plot <- model_data_kec_month %>% 
  group_by(kecid) %>% 
  group_by(treatgroup = factor(control_dummy_median_nvms), yearmonth) %>% 
  filter(year < 2009) %>% 
  summarise(mean_aid = mean(number_of_new_projects, na.rm = T)) %>%
  group_by(yearmonth) %>% 
  ggplot(., aes(x = yearmonth, y = mean_aid, 
                group = treatgroup, 
                fill = treatgroup, 
                color = treatgroup )) +
  geom_point(alpha = 0.4, size = 1.5) +
  geom_line(alpha = 0.4, size = 1) +
  geom_vline(xintercept = as.Date("2006-12-01"), linetype = "dashed") +
  geom_smooth(level = 0.9, se = F, size = 1) +
  geom_label(data = tibble(x = as.Date("2006-12-01"), y = 1.5, label = "Election Month"), 
             aes(x = x, y = y, label = label), inherit.aes = F) + 
  scale_color_brewer("", labels = c("No GAM Area (Control Group)", "GAM Area (Treatment Group)"), palette = "Set1") +
  scale_fill_brewer("", labels = c("No GAM Area (Control Group)", "GAM Area (Treatment Group)"), palette = "Set1") +
  theme_bw() +
  labs(y = "Average Number of New Aid Projects/Month", x = "") +
  theme(panel.grid = element_blank(), 
        legend.position = "bottom")

ggsave(plot = did_plot, filename = "./figures/appendix_figure_6_parallel.pdf", 
       width = 8, height = 5, scale = 0.9)


# Spatial Lag Models ------------------------------------------------------

mod_original  <- felm(log(number_of_new_projects + 1) ~ 
                        pre_election_2006:control_dummy_median_nvms 
                      |  kecid + month | 0 | kecid , 
                      data = model_data_kec_month,
                      keepCX = TRUE  )

aid_mod_sp1 <-  felm(log(number_of_new_projects + 1) ~ 
                       pre_election_2006:control_dummy_median_nvms +
                       log(spatial_lag_tmin1 +1)
                     |  kecid + yearmonth_fct | 0 | kecid , 
                     data = model_data_kec_month   )

mod_sp_50 <- feols(log(number_of_new_projects + 1) ~ 
                     pre_election_2006:control_dummy_median_nvms 
                   |  kecid + month , 
                   data = model_data_kec_month, 
                   vcov = vcov_conley(lat = "latitude", lon = "longitude", 
                                      cutoff = "50"))

mod_sp_100 <- feols(log(number_of_new_projects + 1) ~ 
                      pre_election_2006:control_dummy_median_nvms 
                    |  kecid + month , 
                    data = model_data_kec_month, 
                    vcov = vcov_conley(lat = "latitude", lon = "longitude", 
                                       cutoff = "100"))

mod_sp_200 <- feols(log(number_of_new_projects + 1) ~ 
                      pre_election_2006:control_dummy_median_nvms 
                    |  kecid + month , 
                    data = model_data_kec_month, 
                    vcov = vcov_conley(lat = "latitude", lon = "longitude", 
                                       cutoff = "200"))


## Spatial models table -----

gm <- modelsummary::gof_map %>% 
  filter(grepl("adj\\.r|obs", raw))

rows <- tibble( 
  term = 1:5, 
  "Geographical unit" = rep("Sub-District", 5), 
  "Unit FE" = rep("Yes", 5), 
  "Month FE" = rep("Yes", 5), 
  "Conley SE spatial cluster" = c("-", "-", "50", "100", "200")
) %>% 
  pivot_longer(-term) %>% 
  pivot_wider(names_from = term)

attr(rows, "position") <- c(5, 6, 7, 8, 9)

model_list <- list(mod_original,
                   aid_mod_sp1,
                   mod_sp_50,
                   mod_sp_100,
                   mod_sp_200)

names(model_list) <- paste0("(", as.character(rep(1:5)), ")")

modelsummary(models = model_list, 
             stars = T, 
             coef_map = c("log(spatial_lag_tmin1 + 1)" = "Spatial Lag (t-1)", 
                          "pre_election_2006:control_dummy_median_nvms" = "GAM Area X Pre-2006 Election" ), 
             gof_map = gm, 
             fmt = 2,
             output = "default", 
             add_rows = rows) %>% 
  add_header_above(c(" " = 1, 
                     "No. of New Projcts (log + 1)" = 5)) %>% 
  kableExtra::kable_styling(latex_options = c("hold_position", "scale_down")) %>% 
  kableExtra::save_kable(file = here("./tables/appendix_table_9_spatial.html"))


# Data validity btw ARLS & NVMS -------------------------------------------

data_validity_nvms_arls <- ggplot(model_data_kec_month %>% 
                                    filter(yearmonth == min(yearmonth)),
       aes(x = factor(gam_support_arls),
           y = log(control_nvms + 1))) +
  geom_jitter(size = 2, alpha = 0.6, 
              width = 0.2) +
  geom_boxplot(alpha = 0.6, 
               outlier.alpha = 0) +
  labs(x = "Number of villages in Kecamatan with \nGAM Support == 1 as measured in the village survey", 
       y = "GAM Support as measured by SNPK (log + 1)") +
  theme_bw() + 
  theme(panel.grid = element_blank())

ggsave(plot = data_validity_nvms_arls, 
       filename = "./figures/appendix_figure_5_data_validity.pdf", 
       width = 10, height = 5)


# Survey models for 2007 ----------------------------------------------------------------

baseline_lpm_07 <- felm(election_tsunami_07 ~
                                       gam_support 
                                     | H001a_regency | 0 | H001c_subdistrict, 
                                     data = vh_clean )

covariate_lpm_07 <- felm(election_tsunami_07 ~
                                        gam_support +
                                        # structure
                                        H060i_grp_kpa  +
                                        # demographics
                                        I(tsunami_affected > 0) +
                                        log(H027_vill_pop) +
                                        log(H032a_perc_aceh+1) +
                                        factor(H036_idp) +
                                        H048_poorHH +
                                        H053a_primary98 +
                                        number_of_vh +
                                        I(H110_conf_victims0105 > 0) +
                                        # interviewer FE
                                        factor(H147_lang) +
                                        factor(Hq152_r_6_KPA)
                                      | H001a_regency | 0 | H001c_subdistrict, 
                                      data = vh_clean )

# GLM

baseline_logit_07 <- feglm(election_tsunami_07 ~
                                          gam_support 
                                        | H001a_regency ,
                                        cluster = "H001c_subdistrict",
                                        family = "binomial",
                                        data = vh_clean )

covariate_logit_07 <- feglm(election_tsunami_07 ~
                                           gam_support +
                                           # structure
                                           H060i_grp_kpa  +
                                           # demographics
                                           I(tsunami_affected > 0) +
                                           log(H027_vill_pop) +
                                           log(H032a_perc_aceh+1) +
                                           factor(H036_idp) +
                                           H048_poorHH +
                                           H053a_primary98 +
                                           number_of_vh +
                                           I(H110_conf_victims0105 > 0) +
                                           # interviewer FE
                                           factor(H147_lang) +
                                           factor(Hq152_r_6_KPA)
                                         | H001a_regency ,
                                         cluster = "H001c_subdistrict",
                                         family = "binomial",
                                         data = vh_clean )


## output table ------------------------------------------------------------

gm <- modelsummary::gof_map %>% 
  filter(grepl("adj\\.r|AIC|obs", raw))

rows <- tibble( 
  term = 1:4, 
  # "West Coast Excluded?" = rep("Yes", 8), 
  "Kab. FE" = rep("Yes", 4), 
  "Covariates" =  rep(c("No", "Yes"), 2)) %>% 
  pivot_longer(-term) %>% 
  pivot_wider(names_from = term)

attr(rows, "position") <- c(3, 3, 3)

model_list <- list(baseline_lpm_07, 
                   covariate_lpm_07,
                   baseline_logit_07, 
                   covariate_logit_07)

names(model_list) <- paste0("(", as.character(rep(1:4)), ")")

modelsummary(models = model_list, 
             stars = T, 
             coef_map = c("gam_support" = "Majority in village supported GAM" ), 
             gof_map = gm, 
             output = "default", 
             add_rows = rows) %>% 
  add_header_above(c(" " = 1, 
                     "OLS" = 2, 
                     "Logit" = 2)) %>% 
  add_header_above(c(" " = 1, 
                     "Did village receive tsunami aid in 2007?" = 4)) %>% 
  kableExtra::save_kable(file = here("./tables/appendix_table_11_survey_2007.html"))


# Without West Coast ------------------------------------------------------

# defining West Coast Kabupaten
west_coast_kab <- c("1115", 
                    "1107", 
                    "1116", 
                    "1108")

# 2005

baseline_lpm_05_no_westcoast <- felm(election_tsunami_05 ~
                                       gam_support 
                                     | H001a_regency | 0 | H001c_subdistrict, 
                                     data = vh_clean %>% 
                                       filter(!(H001a_regency %in% west_coast_kab)))

covariate_lpm_05_no_westcoast <- felm(election_tsunami_05 ~
                                        gam_support +
                                        # structure
                                        H060i_grp_kpa  +
                                        # demographics
                                        I(tsunami_affected > 0) +
                                        log(H027_vill_pop) +
                                        log(H032a_perc_aceh+1) +
                                        factor(H036_idp) +
                                        H048_poorHH +
                                        H053a_primary98 +
                                        number_of_vh +
                                        I(H110_conf_victims0105 > 0) +
                                        # interviewer FE
                                        factor(H147_lang) +
                                        factor(Hq152_r_6_KPA)
                                      | H001a_regency | 0 | H001c_subdistrict, 
                                      data = vh_clean %>% 
                                        filter(!(H001a_regency %in% west_coast_kab)))

# GLM

baseline_logit_05_no_westcoast <- feglm(election_tsunami_05 ~
                                          gam_support 
                                        | H001a_regency ,
                                        cluster = "H001c_subdistrict",
                                        data = vh_clean %>% 
                                          filter(!(H001a_regency %in% west_coast_kab)))

covariate_logit_05_no_westcoast <- feglm(election_tsunami_05 ~
                                           gam_support +
                                           # structure
                                           H060i_grp_kpa  +
                                           # demographics
                                           I(tsunami_affected > 0) +
                                           log(H027_vill_pop) +
                                           log(H032a_perc_aceh+1) +
                                           factor(H036_idp) +
                                           H048_poorHH +
                                           H053a_primary98 +
                                           number_of_vh +
                                           I(H110_conf_victims0105 > 0) +
                                           # interviewer FE
                                           factor(H147_lang) +
                                           factor(Hq152_r_6_KPA)
                                         | H001a_regency ,
                                         cluster = "H001c_subdistrict",
                                         family = "binomial",
                                         data = vh_clean %>% 
                                           filter(!(H001a_regency %in% west_coast_kab)))

# 2006

baseline_lpm_06_no_westcoast <- felm(election_tsunami_06 ~
                                       gam_support 
                                     | H001a_regency | 0 | H001c_subdistrict, 
                                     data = vh_clean %>% 
                                       filter(!(H001a_regency %in% west_coast_kab)))

covariate_lpm_06_no_westcoast <- felm(election_tsunami_06 ~
                                        gam_support +
                                        # structure
                                        H060i_grp_kpa  +
                                        # demographics
                                        I(tsunami_affected > 0) +
                                        log(H027_vill_pop) +
                                        log(H032a_perc_aceh+1) +
                                        factor(H036_idp) +
                                        H048_poorHH +
                                        H053a_primary98 +
                                        number_of_vh +
                                        I(H110_conf_victims0105 > 0) +
                                        # interviewer FE
                                        factor(H147_lang) +
                                        factor(Hq152_r_6_KPA)
                                      | H001a_regency | 0 | H001c_subdistrict, 
                                      data = vh_clean %>% 
                                        filter(!(H001a_regency %in% west_coast_kab)))

# GLM

baseline_logit_06_no_westcoast <- feglm(election_tsunami_06 ~
                                          gam_support 
                                        | H001a_regency ,
                                        cluster = "H001c_subdistrict",
                                        family = "binomial",
                                        data = vh_clean %>% 
                                          filter(!(H001a_regency %in% west_coast_kab)))


covariate_logit_06_no_westcoast <- feglm(election_tsunami_06 ~
                                           gam_support +
                                           # structure
                                           H060i_grp_kpa  +
                                           # demographics
                                           I(tsunami_affected > 0) +
                                           log(H027_vill_pop) +
                                           log(H032a_perc_aceh+1) +
                                           factor(H036_idp) +
                                           H048_poorHH +
                                           H053a_primary98 +
                                           number_of_vh +
                                           # factor(H053f_mos98) +
                                           I(H110_conf_victims0105 > 0) +
                                           # interviewer FE
                                           factor(H147_lang) +
                                           factor(Hq152_r_6_KPA)
                                         | H001a_regency ,
                                         cluster = "H001c_subdistrict",
                                         family = "binomial",
                                         data = vh_clean %>% 
                                           filter(!(H001a_regency %in% west_coast_kab)))


## output table ------------------------------------------------------------

gm <- modelsummary::gof_map %>% 
  filter(grepl("adj\\.r|AIC|obs", raw))

rows <- tibble( 
  term = 1:8, 
  "West Coast Excluded?" = rep("Yes", 8), 
  "Kab. FE" = rep("Yes", 8), 
  "Covariates" =  rep(c("No", "Yes"), 4)) %>% 
  pivot_longer(-term) %>% 
  pivot_wider(names_from = term)

attr(rows, "position") <- c(3, 3, 3)

model_list <- list(baseline_lpm_05_no_westcoast, 
                   covariate_lpm_05_no_westcoast,
                   baseline_logit_05_no_westcoast, 
                   covariate_logit_05_no_westcoast, 
                   baseline_lpm_06_no_westcoast,
                   covariate_lpm_06_no_westcoast,
                   baseline_logit_06_no_westcoast,  
                   covariate_logit_06_no_westcoast)

names(model_list) <- paste0("(", as.character(rep(1:8)), ")")

modelsummary(models = model_list, 
             stars = T, 
             coef_map = c("gam_support" = "Majority in village supported GAM" ), 
             gof_map = gm, 
             output = "default", 
             add_rows = rows) %>% 
  add_header_above(c(" " = 1, 
                     "OLS" = 2, 
                     "Logit" = 2, 
                     "OLS" = 2, 
                     "Logit" = 2)) %>% 
  add_header_above(c(" " = 1, 
                     "Did village receive tsunami aid in 2005?" = 4, 
                     "Did village receive tsunami aid in 2006?" = 4)) %>% 
  kableExtra::kable_styling(latex_options = c("hold_position", "scale_down")) %>% 
  kableExtra::save_kable(file = here("./tables/appendix_table_10_westcoast_excluded.html"))


# 2006 Candidate Panel -----------------------------------------------------

candidate_grid <- ggplot(election_2006 %>% 
                           filter(Year == 2006) %>% 
                           mutate(candidates_name = gsub("/", "\n", candidates_name )) %>% 
                           mutate(rebel_candidate = ifelse(grepl("Irwandi", candidates_name), TRUE, FALSE)) %>% 
                           group_by(candidates_name) %>% 
                           # sort candidate names by strength of correlation to control
                           mutate(cor_strength = coef(lm(vote_share ~ control_nvms))[2]) %>% 
                           ungroup() %>% 
                           mutate(candidates_name = forcats::fct_reorder(candidates_name, 
                                                                         cor_strength, .desc = T)), 
                         aes(x = control_nvms, 
                             y = vote_share, 
                             color = rebel_candidate)) +
  geom_point(size = 2.8, alpha = 0.5) + 
  geom_smooth(method = "lm", se = F) +
  facet_wrap(~ candidates_name) +
  scale_color_manual("", 
                     label = c("No rebel candidate", 
                               "Rebel candidate"), 
                     values = c( "darkgrey","firebrick")) +
  theme_bw() +
  labs(x = "GAM support in Kabupaten", y = "Vote share in 2006 elections") +
  theme(panel.grid = element_blank(), 
        legend.position = "none")

ggsave(candidate_grid, filename = here("figures/appendix_figure_3_2006_results.pdf"), 
       height = 7.5, width = 9, scale = 0.8)


# Additional implications survey data: Different aid types -------------------------------

# testing other aid types

bra_kdp_conflictvictim_mod  <- felm(election_bra_kdp_06  ~
                                      factor(H136_support0105 == 1) +
                                      # structure
                                      H060i_grp_kpa  +
                                      # demographics
                                      I(tsunami_affected > 0) +
                                      log(H027_vill_pop) +
                                      log(H032a_perc_aceh+1) +
                                      factor(H036_idp) +
                                      H048_poorHH +
                                      H053a_primary98 +
                                      number_of_vh +
                                      I(H110_conf_victims0105 > 0) +
                                      # interviewer FE
                                      factor(H147_lang) +
                                      factor(Hq152_r_6_KPA)
                                    | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

gam_livelihood_aid_mod <- felm(election_gam_06  ~
                                 factor(H136_support0105 == 1) +
                                 # structure
                                 H060i_grp_kpa  +
                                 # demographics
                                 I(tsunami_affected > 0) +
                                 log(H027_vill_pop) +
                                 log(H032a_perc_aceh+1) +
                                 factor(H036_idp) +
                                 H048_poorHH +
                                 H053a_primary98 +
                                 number_of_vh +
                                 I(H110_conf_victims0105 > 0) +
                                 # interviewer FE
                                 factor(H147_lang) +
                                 factor(Hq152_r_6_KPA)
                               | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

kdp_aid_mod <- felm(election_kdp_06  ~
                      factor(H136_support0105 == 1) +
                      # structure
                      H060i_grp_kpa  +
                      # demographics
                      I(tsunami_affected > 0) +
                      log(H027_vill_pop) +
                      log(H032a_perc_aceh+1) +
                      factor(H036_idp) +
                      H048_poorHH +
                      H053a_primary98 +
                      number_of_vh +
                      I(H110_conf_victims0105 > 0) +
                      # interviewer FE
                      factor(H147_lang) +
                      factor(Hq152_r_6_KPA)
                    | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

tsunami_aid_mod <- felm(election_tsunami_06 ~
                          factor(H136_support0105 == 1) +
                          # structure
                          H060i_grp_kpa  +
                          # demographics
                          I(tsunami_affected > 0) +
                          log(H027_vill_pop) +
                          log(H032a_perc_aceh+1) +
                          factor(H036_idp) +
                          H048_poorHH +
                          H053a_primary98 +
                          number_of_vh +
                          I(H110_conf_victims0105 > 0) +
                          # interviewer FE
                          factor(H147_lang) +
                          factor(Hq152_r_6_KPA)
                        | H001a_regency | 0 | H001c_subdistrict, data = vh_clean)

all_aid_types_model <- bind_rows(broom::tidy(tsunami_aid_mod) %>% mutate(aid_type = "Tsunami Aid"), 
                                 broom::tidy(kdp_aid_mod) %>% mutate(aid_type = "KDP"), 
                                 broom::tidy(gam_livelihood_aid_mod) %>% mutate(aid_type = "GAM Livelihood\nAssistance"), 
                                 broom::tidy(bra_kdp_conflictvictim_mod) %>% mutate(aid_type = "BRA-KDP Conflict Victim"))

all_aid_types_model <- all_aid_types_model %>% 
  filter(grepl("H136_support", term)) %>% 
  mutate(aid_type = forcats::fct_reorder(aid_type, estimate))

aid_types_survey_plot <- ggplot(all_aid_types_model, aes(x = aid_type, y = estimate)) + 
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error, 
                    ymax = estimate + 1.96 * std.error), 
                width = 0) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_bracket(xmin = "BRA-KDP Conflict Victim", xmax =  "KDP", 
               y.position = 0.14, label = "No discretionary control", 
               label.size = 2.9) +
  geom_bracket(xmin = "Tsunami Aid", xmax =  "GAM Livelihood\nAssistance", 
               y.position = 0.14, label = "Discretionary control", 
               label.size = 2.9) +
  labs(y = "Coefficient estimate", 
       x = "Effect of village wartime GAM support on village receiving aid of type ... in 2006") +
  theme_bw() + 
  theme(panel.grid = element_blank()) +
  scale_y_continuous(breaks = seq(-0.2,0.2, 0.05), limits = c(-0.15, 0.18))

ggsave(aid_types_survey_plot, filename = here("figures/appendix_figure_2_alt_test_H3.pdf"), 
       width = 9, height = 6, scale = 0.8)


# Add. impl. RAN data: aid sector -------------------------------------------

# Project Count as DV

mod_sector_spatial_count <- felm(log(sector_spatial_count + 1) ~ 
                                   pre_election_2006:control_dummy_median_nvms 
                                 |  kecid +yearmonth_fct | 0 | kecid , 
                                 data = model_data_kec_month )

mod_sector_unall_count <- felm(log(sector_unall_count + 1) ~ 
                                 pre_election_2006:control_dummy_median_nvms
                               |  kecid +yearmonth_fct | 0 | kecid , 
                               data = model_data_kec_month )

mod_sector_social_count <- felm(log(sector_social_count + 1) ~ 
                                  pre_election_2006:control_dummy_median_nvms
                                |  kecid +yearmonth_fct | 0 | kecid , 
                                data = model_data_kec_month )

mod_sector_admin_count <- felm(log(sector_admin_count + 1) ~ 
                                 pre_election_2006:control_dummy_median_nvms
                               |  kecid +yearmonth_fct | 0 | kecid , 
                               data = model_data_kec_month )

mod_sector_econ_count <- felm(log(sector_econ_count + 1) ~ 
                                pre_election_2006:control_dummy_median_nvms
                              |  kecid +yearmonth_fct | 0 | kecid , 
                              data = model_data_kec_month )

mod_sector_educ_count <- felm(log(sector_educ_count + 1) ~ 
                                pre_election_2006:control_dummy_median_nvms
                              |  kecid +yearmonth_fct | 0 | kecid , 
                              data = model_data_kec_month )

mod_sector_health_count <- felm(log(sector_health_count + 1) ~ 
                                  pre_election_2006:control_dummy_median_nvms
                                |  kecid +yearmonth_fct | 0 | kecid , 
                                data = model_data_kec_month )

mod_sector_infra_count <-  felm(log(sector_infra_count + 1) ~ 
                                  pre_election_2006:control_dummy_median_nvms
                                |  kecid +yearmonth_fct | 0 | kecid , 
                                data = model_data_kec_month )

mod_sector_institution_count <- felm(log(sector_institution_count + 1) ~ 
                                       pre_election_2006:control_dummy_median_nvms
                                     |  kecid +yearmonth_fct | 0 | kecid , 
                                     data = model_data_kec_month )

mod_sector_religion_count <- felm(log(sector_religion_count + 1) ~ 
                                    pre_election_2006:control_dummy_median_nvms
                                  |  kecid +yearmonth_fct | 0 | kecid , 
                                  data = model_data_kec_month )

models_sector_count <- bind_rows(broom::tidy(mod_sector_spatial_count) %>% mutate(model = "Spatial Planning & Environment"),
                                 broom::tidy(mod_sector_unall_count) %>% mutate(model = "Unallocated"),
                                 broom::tidy(mod_sector_social_count) %>% mutate(model = "Social"),
                                 broom::tidy(mod_sector_admin_count) %>% mutate(model = "Administrative & Other"),
                                 broom::tidy(mod_sector_econ_count) %>% mutate(model = "Economic Development"),
                                 broom::tidy(mod_sector_educ_count) %>% mutate(model = "Education"),
                                 broom::tidy(mod_sector_health_count) %>% mutate(model = "Health"),
                                 broom::tidy(mod_sector_infra_count) %>% mutate(model = "Infrastructure"),
                                 broom::tidy(mod_sector_institution_count) %>% mutate(model = "Institutional Development"),
                                 broom::tidy(mod_sector_religion_count) %>% mutate(model = "Religion"))

models_sector_count <- models_sector_count %>% 
  mutate(model = str_wrap(model, 10)) %>% 
  mutate(model = forcats::fct_reorder(model, estimate * -1))

coefplot_sector_count <- ggplot(models_sector_count, aes(x = model, 
                                                         y = estimate)) +
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error, 
                    ymax = estimate + 1.96 * std.error), 
                width = 0) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Sector", y = "Coefficient Estimate") +
  theme_bw() + 
  theme(panel.grid = element_blank()) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

ggsave(coefplot_sector_count, filename = here("./figures/appendix_figure_1_alt_test_H2.pdf"), 
       height = 4, 
       width = 11, 
       scale = 0.8)
